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ABSTRACT 

We present a test of the emission statistics of active galactic nuclei (AGN), probing the connection 
between red-noise temporal power spectra and multi-modal flux distributions known from observations. 
We simulate AGN lightcurves under the assumption of uniform stochastic emission processes for different 
power-law indices of their power spectra. For sufficiently shallow slopes (power-law indices ft 1), the 
flux distributions (histograms) of the resulting lightcurves are approximately Gaussian. For indices 
corresponding to steeper slopes (/3 J> 1), the flux distributions become multi- modal. This finding 
disagrees systematically with results of recent mm/radio observations. Accordingly, we conclude that 
the emission from AGN does not necessarily originate from uniform stochastic processes even if their 
power spectra suggest this. Possible mechanisms are transitions between different activity states and/or 
the presence of multiple, spatially disconnected, emission regions. 

Key words : galaxies: active — radiation mechanisms: general — methods: statistical 



I. INTRODUCTION 

There is general consensus that active galactic nuclei 
(AGN) are powered by accretion of matter onto super- 
massive (M. « 1O 6 - 1O M ) black holes located in the 
centers of most, if not all, galaxies (see e.g. Peterson 
1997, Krolik 1999, or Beckmann & Shrader 2012, for 
reviews). AGN are known for strong temporal variabil- 
ity of their luminosity. Especially at radio frequencies, 
variability appears to be related to shocks propagating 
through jets (e.g. Marscher & Gear 1985). At optical, 
ultraviolet, and higher frequencies changes in accretion 
flows and accretion disk structure are supposed to fur- 
ther modulate the luminosity (e.g. De Vries et al. 2005; 
Czerny 2006). 

To understand the underlying physics of the variabil- 
ity, several studies aimed at identifying characteristic 
variability timescales. Timescales probed by observa- 
tions range from few kiloseconds (Benlloch et al. 2001; 
using X-ray observations, Schbdel et al. 2007; using ra- 
dio observations) to tens of years (Hovatta et al. 2007, 
2008; using radio observations). Other studies ana- 
lyzed the temporal power spectra - the square mod- 
uli of the Fourier transforms - of AGN lightcurves. 
Press (1978) pointed out that power spectra of AGN 
lightcurves follow power laws Af<xf~^ with j3 > 0, i.e. 
red noise laws; here Af denotes the power spectral am- 
plitude for a given sampling frequency /.* Cases com- 
monly encountered in statistics are white noise (/3 = 0) , 
random walk noise (0 = 2), and (intermediate) flicker 
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noise (f3 = 1). In many cases observed, AGN power 
spectra can be described by 1// flicker noise (Lawrence 
et al. 1987). Red- noise power spectra are analyzed 
in depth by Lawrence & Papadikas (1993) who used 
12 high-quality "long look" X-ray lightcurves of AGN 
from the EXOSAT data base, as well as by Uttley et al. 
(2002) who found observational evidence for a flatten- 
ing of AGN X-ray power spectra at timescales longer 
than few months. At optical wavelengths, De Vries 
et al. (2005) found evidence for the variability of quasar 
emission following red-noise laws on timescales as long 
as approximately 40 years. 

Whether the emission processes of AGN are indeed 
uniform stochastic has been debated hotly especially 
for the case of Sagittarius A*, the supermassive black 
hole located at the center of the Milky Way. Some 
studies argued that power spectra of its near-infrared 
lightcurves are consistent with a single red noise law 
(Meyer et al. 2008; Do et al. 2009). Other studies found 
indication that there are multiple states of emission, 
possibly including a quasi-periodic modulation of the 
source flux (Gillessen et al. 2006, Trippe et al. 2007, 
Dodds-Eden et al. 2011). 

Recent millimeter/radio observations (Trippe et al. 
2011) indicate that the emission statistics of AGN is 
more complex than generally assumed. Whereas the 
lightcurves of the sample of Trippe et al. (2011) showed 
the known red noise power spectra, the flux distribu- 
tions (histograms) showed clear bi- or multi-modality. 
This result was somewhat unexpected as power spectra 
with (3 < 1, as observed by Trippe et al. (2011), can be 
understood as being generated by lightcurves with flux 
values drawn from uni-modal - approximately Gaus- 
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sian - distributions. The observations were interpreted 
qualitatively as indication for the presence of distinct 
emission states. 

A quantitative test of this interpretation requires a 
careful evaluation of the relation between the slope of 
the power spectrum and the corresponding flux distri- 
bution. Red noise power spectra imply higher variabil- 
ity at longer timescales, and realistic lightcurves as well 
as the corresponding power spectra are noisy. Accord- 
ingly, extreme flux values or multi-modal flux distribu- 
tions may occur under certain conditions even if not 
expected in general. 

II. ANALYSIS 

A. Simulated Lightcurves 

Monte Carlo simulations provide a convenient way to 
probe the relations between red noise power spectra 
and the underlying flux distributions. We generate ar- 
tificial Fourier transforms of lightcurves with (3 = 0, 
0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, and 2 using the al- 
gorithm introduced by Timmer & Kdnig (1995). For 
each sampling frequency /, we draw two random num- 
bers from Gaussian distributions for the real part and 
the imaginary part, respectively. Power-law noise with 
slope —f3 is generated by multiplying both values with 
^dividual values are stored in an array 
of complex numbers; this array corresponds to the 
complex Fourier transform of our artificial lightcurve. 
We calculate the normalized temporal power spectrum, 
also known as periodogram, by taking the square of the 
absolute value of the array. A periodogram is given by 



|FH| 2 



( s(t) cos ^ + ( X/ s ® sin wt ) 



F(—fi) = F*(fi)i where the operator * denotes com- 
plex conjugation (Timmer & Kdnig 1995). For each 
value of (3 we simulate TV' - meaning the number of 
trials - lightcurves. 

Each artificial lightcurve is composed of 20 000 data 
points initially. However, in order to circumvent intrin- 
sic symmetries of the Fourier transform process, it is 
necessary to restrict the analysis to a subset of - ad- 
jacent - data points taken from one half of the time 
scries. A priori, the choice of a certain number of data 
points should not influence the properties of the result- 
ing flux distributions; this is due to the use of power- 
law noise and the fact that power-laws are scale free. 
Nevertheless, to minimize the risk of selection effects, 
we perform two separate analysis runs: 

> Run A uses N = 9 000 data points per time series, 
the number of trials is N' = 100. 

> Run B uses N = 2 000 data points per time series, 
the number of trials is N' = 500. 

In order to ease comparison among different realiza- 
tions of lightcurves, we use the normalization 



S(U) - S(t) 



(3) 



(i) 



where S(t) is the lightcurve, S(t) denotes the mean 
value of S(t), and a is the standard deviation of S(t). 
Accordingly, our time series are normalized to zero 
mean and unity standard deviation. 

We check that our choices of N do not alter the 
values of (3 - which were applied to time series with 
20 000 data points initially - systematically. For this 
check, we obtain the actual (3 values by calculating for 
each lightcurve the normalized Scargle periodogram 



(Schuster 1898) where S(t) is the signal in time do- 
main, F(w) is the Fourier transform, and lo = 2nf. 
By construction, a periodogram of our simulated sig- 
nal follows the relation Ajocf~P. The results agree 
with power spectra of actual observations (see Fig. 1). 

We obtain simulated lightcurves S(t) by taking the 
inverse Fourier transform of our complex arrays. Here 
we use the inverse Fourier transform of the FFT (Fast 
Fourier Transform) procedure in IDL.t The inverse 
Fourier transform is given by 



N/2 

S W = E F(f)eM^ft/N] 

f=-N/2 



(2) 



where t is the time, / denotes the sampling frequency, 
F(f) are the values in frequency domain and N is the 
total number of array elements. To obtain a real valued 
time series, we construct our complex arrays such that 



J2i cos 2 2tt fti 



J^i sin 2 : 2tt fU 



(4) 



(Scargle 1982). Here Af is the amplitude of the peri- 
odogram evaluated at frequency /, which is equivalent 
to the power spectral density;* Si are the flux value 
and ti the time of the ith data point, respectively; and 
a 2 is the variance of the data. The base frequency is 
/ m ; n ; the ith frequency is i times of / m in, i-e. / = 

/minj 2/ m i n , 3 /mini • • • ; /max with /min — /max — 

N/ (2T). Here T is the total observation time and N is 
the number of data points (9 000 for run A and 2 000 
for run B). We confirm that input values and actual 
values of (3 agree within w 2%. 



tFor the analysis of observational data, the Scargle periodogram 
is preferable over the standard periodogram as it can be applied 
to data with arbitrary sampling and has a well-understood sta- 
tistical behavior. 



t Interactive Data Language, ITT Exelis Inc., McLean (Virginia). 
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Fig. 1. — Simulated power spectra with /3 = 0, 1,2 (from left to right) in logarithmic scales. Units are arbitrary. These 
artificial power spectra agree with actual observations (compare, e.g. Timmer & Konig 1995, or Trippe et al. 2011). 



B. Exploring the Flux Distributions 

For each artificial lightcurve, we obtain its flux distri- 
bution by calculating the histogram of the flux values; 
example results are provided in Fig. 2. By construc- 
tion, these flux distributions correspond to the results 
of single, uniform stochastic emission processes. De- 
pending on the value of /3, we find uni-modal, approx- 
imately Gaussian, as well as multi-modal profiles. For 
(3 <; 1.25, the flux distributions agree qualitatively with 
noisy Gaussian distributions. For (3 J> 1.25, the distri- 
butions become flattened, and multi-modality becomes 
evident with increasing (3. 

In order to check for deviations from uni-modality 
quantitatively, we probe the agreement of flux his- 
tograms with Gaussian distributions by means of a x 2 
test (see also Fig. 2) with 



N" 

E 



Ern 2 



(5) 



(e.g. Andrae et al. 2010) where Si denotes the ith value 
of the given distribution, Ei denotes the value expected 
theoretically, and Erri is the statistical error of the 
given value. Errors are binomial. For our theoretical 
Gaussian function we use the parametrization 



f(x;n,a 2 ) = Aexp 



(6) 



where A is the amplitude, a 2 is the variance, and \x 
is the mean of the distribution. Because of our nor- 
malization of the simulated flux distributions (Eq. 3), 
we have a 2 — 1 and fi = 0, leaving A as the only free 
parameter. We also obtain the reduced x 2 given by 



Xrcd 



X 



DOF 



(7) 



where DOF is the number of degrees of freedom of the 
model given by 



DOF = TV" - P 



(8) 



(e.g. Andrae et al. 2010) where N" is the number of 
bins of the histogram that have non-zero values and 
P is the number of parameters of the model. In our 
case, P = 1 because the amplitude of the Gaussian 
function is our only model parameter (see Eq. 6). If 
the observed distribution corresponds to the expected 
one within errors, Xrcd ~ 1- We apply x 2 tests to each 
simulated flux distribution for each value of (3; example 
results of these tests are shown in Fig. 2. We observe a 
systematic increase of x 2 0( j w ith increasing /?, suggest- 
ing substantial deviation from Gaussian distributions 
occurring for /? > 1. 

Whereas a x 2 test provides information on the pres- 
ence of deviations between model and data, it does not 
tell us how model and data deviate. We are here specif- 
ically concerned about a "breakup" of a single flux dis- 
tribution into multiple, potentially overlapping, distri- 
butions with increasing values of (3 as illustrated in Fig. 
2. Accordingly, we find distributions that can deviate 
systematically from Gaussians in terms of symmetry as 
well as peakedness (concentration). 

The degree of asymmetry of a probability distribu- 
tion is quantified by the skewness, or third standardized 
moment, 



.91 = 3 



(E[(X-^]f/2 



(9) 



(e.g. Crawley 2005; Press et al. 2007). Here E is the 
expectation value operator, X is the data set, /13 is the 
third moment about the mean /x, and er is the standard 
deviation. If the probability distribution is Gaussian, 
the skewness is zero by construction. For each value 
of j3, we calculate g\ for all N' realizations of flux dis- 
tributions; we present the corresponding histograms in 
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Fig. 2. — Realizations of simulated flux distributions, i.e. histograms of simulated lightcurves (continuous lines), 
together with the best- fitting Gaussian distributions (dashed lines), for all eight non-zero values of j3. Lightcurves, and thus 
their histograms, are normalized to zero mean and unity standard deviation; accordingly, abscissae are in arbitrary count 
units. Ordinates denote the number of data points in a given flux bin. Top group (above the horizontal line): Results for 
RUN A; each histogram is composed of 9 000 flux values. Bottom group (below the horizontal line): Results for run B; each 
histogram is composed of 2 000 flux values. Please note the different axis scales. In each diagram, we quote the \ 2 an( A Xrcd 2 
values for the Gaussian fits. For /3 ;> 1.25, the flux distributions break up into several overlapping distributions. 
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Fig. 3. — Histograms of the skewness of the simulated flux distributions for all eight non-zero values of /3. (Examples for 
flux distributions arc shown in Fig. 2.) Abscissae indicate the (binned) skewness values observed; ordinates show the number 
of simulated flux distributions with a skewness located in a given bin. Top group (above the horizontal line): Results for 
RUN A; each histogram includes the results from 100 trials. Bottom group (below the horizontal line): Results for run B; 
each histogram includes the results from 500 trials. Please note the different axis scales. The average values of skewness 
do not significantly - in terms of standard error of mean - deviate from zero, meaning that the given functions agree with 
Gaussian profiles with respect to symmetry. 
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Fig. 4. — Histograms of the KURTOSIS of the simulated flux distributions for all eight non-zero values of /?. Abscissae 
indicate the (binned) kurtosis values observed; ordinates show the number of simulated flux distributions with a kurtosis 
located in a given bin. Top group (above the horizontal line): Results for run A; each histogram includes the results from 
100 trials. Bottom group (below the horizontal line): Results for run B; each histogram includes the results from 500 trials. 
Please note the different axis scales. The average values are significantly - in terms of standard error of mean - smaller 
than zero for > 1.0, meaning that the observed distributions are too flat compared to Gaussian profiles (platykurtic 
distributions). This signature is expected for multi-modal flux distributions. 
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Table 1. 

Skewness statistics for runs A and B. Deviations from zero are in units of a-g, i.e. the standard error of mean. 



0.0 



0.25 



0.5 



0.75 



1.0 



1.25 



1.5 



1.75 



2.0 



Run A: 
mean 

standard deviation 
standard error 
deviation from (a^) 

Run B: 
mean 

standard deviation 
standard error 
deviation from (cr^) 



0.003 0.003 

0.023 0.023 

0.002 0.002 

1.261 1.335 



0.004 0.005 

0.056 0.056 

0.003 0.003 

1.597 1.996 



0.003 
0.027 
0.003 
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0.003 
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0.047 
0.005 
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0.004 
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0.015 
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0.011 
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0.158 
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0.036 
0.219 
0.022 
1.635 



0.002 
0.258 
0.012 
0.173 



0.055 
0.309 
0.031 
1.774 



0.002 
0.350 
0.016 
0.128 



0.061 
0.373 
0.037 
1.628 



0.055 
0.422 
0.042 
1.302 



0.000 -0.004 

0.421 0.474 

0.019 0.021 

0.000 0.189 



Fig. 3. Within errors, the average values of skewness 
agree with zero for all (3. We note, however, that the 
width of the skewness distributions - quantified by the 
standard deviations of the histograms - increases with 
increasing /?, indicating a stronger scatter of skewness 
values. 

The concentration, or peakedness, of a probability 
distribution is quantified by the kurtosis, or fourth 
standardized moment, 



= m 4 E[(X-rf] 
92 rnj {E[(X - 



(10) 



(e.g. Crawley 2005; Press et al. 2007). Here m 4 de- 
notes the fourth sample moment about the mean, and 
m 2 is the second sample moment about the mean. If 
the probability distribution is Gaussian, the kurtosis 
is zero by construction. For each value of /3, we cal- 
culate §2 for all N 1 realizations of flux distributions; 
we present the corresponding histograms in Fig. 4. 
For (3 < 1 , the average values of kurtosis are in agree- 
ment with zero (within errors), indicating agreement 
with Gaussian flux distributions. For (3 J> 1, the mean 
values of kurtosis deviate significantly from zero; the 
mean values found in runs A and B are in agree- 
ment within (standard) errors. Notably, the average 
kurtosis becomes negative, meaning the simulated flux 
distributions are, in average, platykurtic - they are too 
flat compared to Gaussian profiles. This is indeed the 
signature we expect if a single, uni-modal distribution 
breaks up into multiple overlapping distributions. 

In order to provide a more quantitative assessment, 
we check how much the averages of skewness and kur- 
tosis deviate from zero in units of the standard error of 
the mean, which is given by 



a 



0~x 



(11) 



(Diekhoff 1992). Here N' is - again - the number of 
trials, being 100 for run A and 500 for run B. The 
results are shown in Tables 1 and 2. Table 1 shows that 
the average values of skewness never deviate by more 
than w 2 er— from zero. In contrast, Table 2 shows that 
the average values of kurtosis deviate significantly - 
meaning by more than 5 a— - from zero for j3 ;> 1.0. 

III. DISCUSSION 

Our analysis provides a valuable test of the statisti- 
cal emission properties of active galactic nuclei. We 
explore the relation between the shapes of the flux dis- 
tributions (histograms) of red-noise lightcurves on the 
one hand, and the power-law indices of the temporal 
power spectra (periodograms) of those lightcurves on 
the other hand. For values (3 <J 1, the flux distribu- 
tions are in agreement (see Tables 1, 2) with being 
Gaussian in average. This picture changes gradually 
for (3 ;> 1 : with increasing values of {3 the flux distribu- 
tions tend to break up into multiple, overlapping sub- 
distributions (illustrated in Fig. 2), leading to statisti- 
cally significant deviations from Gaussian distributions 
(Fig. 4; Tables 1, 2). Turning this finding around, we 
may conclude: if the power-law index is smaller than 
unity and observations find multi-modal flux distribu- 
tion then the emission does not originate from a single, 
uniform stochastic emission process. 

Recent observations point out the necessity to an- 
alyze flux distributions carefully when addressing the 
emission statistics of AGN. Dodds-Edcn et al. (2011) 
gave observational evidence for multiple emission states 
of Sgr A*. They concluded that its near-infrared 
flux distribution deviates systematically from realistic 
analytical profiles, especially log-normal distributions. 
This discrepancy can be resolved by assuming a super- 
position of two temporally distinct states of activity: a 
"quiescent" state at low energies, and a "flare" state at 
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Table 2. 

KuRTOSlS statistics for runs A AND B. Deviations from zero are in units of a—, i.e. the standard error of mean. 



(3 0.0 0.25 0.5 0.75 1.0 1.25 1.5 1.75 2.0 

Run A: 
mean 

standard deviation 
standard error 
deviation from (a^) 

Run B: 
mean 

standard deviation 
standard error 
deviation from (ax) 



-0.004 
0.052 
0.005 
0.710 



-0.003 
0.052 
0.005 
0.579 



-0.002 
0.053 
0.005 
0.375 



-0.009 
0.072 
0.007 
1.222 



-0.059 
0.154 
0.015 
3.833 



-0.184 
0.301 
0.030 
6.097 



-0.345 
0.417 
0.042 
8.272 



-0.494 
0.496 
0.050 
9.960 



-0.612 
0.560 
0.056 
10.92 



0.003 
0.114 
0.005 
0.588 



0.001 
0.114 
0.005 
0.196 



-0.003 
0.114 
0.005 
0.588 



-0.021 
0.131 
0.006 
3.585 



-0.085 
0.203 
0.009 
9.363 



-0.217 
0.323 
0.014 
15.02 



-0.381 
0.435 
0.019 
19.59 



-0.525 
0.518 
0.023 
22.66 



-0.635 
0.580 
0.026 
24.48 



high energies. 

Trippc et al. (2011) analyzed the high-frequency ra- 
dio lightcurves of six AGN. Five of the six sources 
showed an unexpected pattern: whereas their peri- 
odograms followed red noise laws (with 0.4 <;/?<; 0.7)§, 
their flux distributions showed strong bi- or even multi- 
modality. This behavior led Trippe et al. (2011) to sus- 
pect that the emission observed originates from distinct 
emission states. One possible explanation is provided 
by temporally distinct states of activity, like quiescent 
and, potentially multiple, flare states (cf. Dodds-Eden 
et al. 2011). Another explanation is the superposi- 
tion of flux from several, causally disconnected emis- 
sion regions - like different shock zones in radio jets. 
Trippe et al. (2011) argued qualitatively that for red- 
noise power spectra with /3 <; 1, the observed flux dis- 
tributions should be uni-modal if the emission origi- 
nates from the same stochastic process. However, this 
argument was never tested quantitatively; such a test 
is now provided by the present study. 

Our statistical analysis indicates that, at least for 
moderate red- noise power- law indices of [3 <J 1, the un- 
derlying time series should be drawn from a uni-modal, 
approximately Gaussian, distribution. Evidently, the 
recent observations by Trippe et al. (2011) suggest that 
the actual picture is more complicated: even if the 
power spectra apparently correspond to uniform red- 
noise processes, the emission can arise from multiple 
components - where these "components" may be dis- 
tinct in time (temporally distinct states of activity) as 
well as in space (like different emission zones in jets). 
A combined analysis of temporal power spectra and 
flux distributions may be a valuable tool for studying 
AGN structures in general: as the analysis is based 
on lightcurves only, it probes the structure of sources 



§When taking into account the measurement noise, the intrin- 
sic values of /3 are slightly larger but still well below unity (see 
Sect. 5.3 of Trippe et al. 2011 for details). 



without actually resolving them spatially. 

Mathematically, red-noise power spectra are due to 
correlations between adjacent data points in the under- 
lying time series (Press 1978). The actual power- law 
index is given by the strength of correlation, ranging 
from no correlation at all - white noise, {3 = - up 
to correlation over infinite times - random walk noise, 
(3 = 2. However, a physical explanation for the oc- 
currence of red-noise time series in many physical and 
astrophysical processes is not yet known, and is subject 
to ongoing studies (Press 1978; Kaulakys & Alaburda 
2009; Kelly, Sobolewska & Siemiginowska 2011). Ac- 
cordingly, new insights and surprises regarding the sta- 
tistical properties of AGN emission should be expected 
for the future. 

IV. CONCLUSIONS 

This study presents a test of the statistical properties 
of emission from active galactic nuclei. Using the algo- 
rithm introduced by Timmer & Konig (1995), we sim- 
ulate lightcurves based on the assumption of uniform 
stochastic red-noise (0 < (3 < 2) emission processes. 
We compare the resulting flux distributions to those of 
recent mm/radio observations by Trippe et al. (2011). 
Our work arrives at the following conclusions: 

1. For power-law indices (3 <; 1, the flux distributions 
(histograms) of our simulated lightcurves agree 
with noisy Gaussian profiles. For (3 ;> 1, the 
flux distributions remain symmetric (zero skew- 
ncss) but become increasingly platykurtic (nega- 
tive kurtosis) with increasing values of ft. In those 
cases, the flux distributions break up into several 
overlapping sub-distributions, making the overall 
distributions multi-modal. 

2. From the behavior of the flux distributions as 
function of (3, we also conclude: if the power-law 
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index is smaller than unity and observations un- 
veil multi-modal flux distributions, then the emis- 
sion does not arise from a single, uniform stochas- 
tic emission process. This relation may be used 
to probe the structure of AGN without resolving 
the target sources spatially. 

3. Comparing our statistical analysis with the recent 
observation of multi- modal mm/radio flux distri- 
butions of five AGN by Trippe et al. (2011), it 
becomes evident that Trippe et al. (2011) indeed 
observed emission from multiple emission states 
- temporally distinct states of activity and/or 
spatially disconnected emission regions - in their 
sources. In general, it appears that AGN emis- 
sion can originate from multiple emission pro- 
cesses even if their simple red-noise power spectra 
suggest otherwise. 

A combined analysis of temporal power spectra and 
flux distributions of AGN can, a priori, be applied to a 
large variety of data sets available, like e.g. the radio 
flux monitoring by Hovatta et al. (2007, 2008). Even 
though a priori simple, it appears that the power of 
such an analysis has largely been overlooked so far. 
Accordingly, we expect new insights from extended sta- 
tistical analyses of AGN lightcurves in the future. 
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